      !< @file
      !> @brief Subrotina para o cálculo do erro da solução do problema inverso.
      !>
      !> Calcula o erro da imagem obtida com a inversão em relação à imagem do modelo de velocidades usado para geração dos sismogramas.

      subroutine calcErro(nx,nz,npml,unit)

      implicit none

      integer :: nx !< número de nós do grid na direção x
      integer :: nz !< número de nós do grid na direção z      
      integer :: npml !< número de nós em cada camada de absorção
      integer :: unit !< unidade onde serão impressos os erros do modelo aproximado em relação ao modelo real

      double precision :: erromedio !< Média dos valores absolutos do erro em cada nó
      double precision :: erromax !< Valor máximo absoluto para o erro
      double precision :: erromediorelat !< Média dos valores absolutos do erro em cada nó dividido pelo valor do modelo exato
      double precision :: erromaxrelat !< Valor máximo absoluto para o erro dividido pelo respectivo valor do modelo exato

      real, dimension(:,:), allocatable :: Modelo !< modelo usado para geração dos sismogramas
      real, dimension(:,:), allocatable :: Minv  !< modelo obtido com a inversão 
      integer :: ierr !< identificação de erros
      integer :: i !< contador na direção x
      integer :: j !< contador na direção z

      integer :: contnohs !< conta os nós usados para cálculo de erro




      allocate(Modelo(nx,nz),stat=ierr)
      allocate(Minv(nx,nz),stat=ierr)


      open(1,file='veloc.bin',status="old",form="unformatted",access="direct",recl=4*Nx*Nz,iostat=ierr)
      if(ierr.eq.0)then
	  read(1,rec=1)((Modelo(i,j),j=1,Nz),i=1,Nx)
	  close(1)
      else
	  write(*,*)'Nao foi possível ler o arquivo veloc.bin :  IOSTAT=',ierr
	  erromedio=-1.0
	  return
      endif

      open(1,file='veloX.bin',status="old",form="unformatted",access="direct",recl=4*Nx*Nz,iostat=ierr)
      if(ierr.eq.0)then
	  read(1,rec=1)((Minv(i,j),j=1,Nz),i=1,Nx)
	  close(1)
      else
	  write(*,*)'Nao foi possível ler o arquivo veloX.bin :  IOSTAT=',ierr
	  erromedio=-1.0
	  return
      endif


      erromedio = 0.0d0
      erromediorelat = 0.0d0
      contnohs = 0
      do i=Npml+1,Nx-Npml
	  do j=1,Nz-Npml
	      erromedio = erromedio + abs(Modelo(i,j)-Minv(i,j))
	      erromediorelat = erromediorelat + abs(Modelo(i,j)-Minv(i,j))/Modelo(i,j)
	      contnohs = contnohs+1
	  enddo
      enddo
      erromedio = erromedio/dfloat(contnohs)
      erromediorelat = erromediorelat/dfloat(contnohs)


     

      erromaxrelat = 0.0d0
      do i=Npml+1,Nx-Npml
	  do j=1,Nz-Npml
	      if(abs(Modelo(i,j)-Minv(i,j))/Modelo(i,j).gt.erromaxrelat)then
		    erromaxrelat = abs(Modelo(i,j)-Minv(i,j))/Modelo(i,j)
		    erromax = abs(Modelo(i,j)-Minv(i,j))
	      endif
	  enddo
      enddo


      write(*,*)'Erro medio: ', erromedio,'m/s  =  ',erromediorelat*100.0,'%'
      write(*,*)'Erro maximo relativo: ', erromaxrelat*100.0d0,'%  =', erromax,'m/s'


      write(unit,*)'#Erro medio: ', erromedio,'m/s  =  ',erromediorelat*100.0,'%'
      write(unit,*)'#Erro maximo relativo: ', erromaxrelat*100.0d0,'%  =', erromax,'m/s'







      deallocate(Modelo)
      deallocate(Minv)




      return
      end subroutine calcErro
    